## replication file for Cook-Fortunato ##
## load libraries and data ##

	library(doBy)
	library(geofacet)
	library(ggplot2)
	library(ggpubr)
	library(haven)
	library(lme4)
	library(stargazer)
	library(tidyr)

	old <- read_dta('ucrOld.dta')
	new <- read_dta('ucrNew.dta')

## variables: ##
	## ucr = outcome, indicates compliance with ucr as defined in main tex
	## faOrd = central covariate, factor decomposed scale of legislative salary, session length, and spending
	## termEnact = central covariate, indicates term limits have been adopted
	## election = indicates agency run by elected sheriff
	## conSheriff = indicates sheriff's office is mentioned in state constitution
	## type = indicates agency organized at city, county, or state level
	## employees = total civilian and officer emplyoee count for the agency-year
	## population = total population governed by agency as identified by FBI
		## note 0 counts are common for state outposts
	## demGov = indicates democratic governor (relative to independent)
	## repGov = indicates republican governor (relative to independent)
	## demLeg = indicates democratic legislative control (relative to divided)
	## repLeg = indicates republican legislative control (relative to divided)
	## black = percentage of state population identified as balck in most recent census
	## state = identifies state housing the focal agency
	## year = identifies year for which ucr request was made
	## ori = 'originating agency identification' agency's unique federal id

## First, study 1, ucr compliance ##
## Note that these data have imputed capacity for years preceding 1974 ##
## and impute missingness on employee counts and population (hence the negatives). ##
## We do not calculate Rubin's errors, however, because there is ##
## effectively 0 variance in key covariate parameter and standard error estimates. ##
## If you are using these data for new research, however, please use the files ##
## labeled "raw data" and re-impute with the associated script. ##

## Table 1 models ##
## 1960-1994 ##
## capacity: reduced form model ##
	squireOld <- lm(ucr ~ faOrd + 
		as.factor(state) + as.factor(year), data = old)
	summary(squireOld)

## term limits: reduced form model ##
	termOld <- lm(ucr ~ termEnact + 
		as.factor(state) + as.factor(year), data = old)
	summary(termOld)

## all covariates ##
	fullOld <- lm(ucr ~ (faOrd + termEnact) * election + 
		type + conSheriff + scale(employees) + scale(population) + 
		termEnact + demGov + repGov + black + 
		as.factor(state) + as.factor(year)*(demLeg + repLeg), data = old)
	summary(fullOld)

## 1995-2017 ##
## capacity: reduced form model ##
	squireNew <- lm(ucr ~ faOrd +
		as.factor(state) + as.factor(year), data = new)
	summary(squireNew)

## term limits: reduced form model ##
	termNew <- lm(ucr ~ termEnact +
		as.factor(state) + as.factor(year), data = new)
	summary(termNew)

## all covariates ##
	fullNew <- lm(ucr ~ (faOrd + termEnact) * election + 
		type + conSheriff + scale(employees) + scale(population) + 
		termEnact + demGov + repGov + black + 
		as.factor(state) + as.factor(year)*(demLeg + repLeg), data = new)
	summary(fullNew)

## table 1 output ##
	stargazer(squireOld, termOld, fullOld, squireNew, termNew, fullNew)

## now study 2, underporting analyses ##
	data <- read_dta('underReport.dta')

## variables: ##
	## police_homicide = outcome, count of killings by police in ucr
	## police_homicide_wp = outcome, count of killings by police by washington post
	## police_homicide_guardian = outcome, count of killings by police by guardian
	## police_homicide_mpv = outcome, count of killings by police by mapping police violence
	## police_homicide_kbp = outcome, count of killings by police by killed by police
	## police_homicide_fatal = outcome, count of killings by police by fatal encounters
	## std_squire_capacity = central covariate, same factor decomposed scale above, normalized
	## termEnact = central covariate, indicates term limits have been adopted
	## invLaw = central covariate, indicates police homocide investigation policy active

	## std_med_income = state year median income, normalized
	## std_employees = state-year avergage agency civilian and officer employee counts, normalized
	## std_police_ideal = average ideal point (bonica data) from all contributors in policing, normalized
	## demGov = indicates democratic governor (relative to independent)
	## repGov = indicates republican governor (relative to independent)
	## demLeg = indicates democratic legislative control (relative to divided)
	## repLeg = indicates republican legislative control (relative to divided)
	## std_black = percentage of state population identified as balck in most recent census, normalized
	## state = identifies state housing the focal agency
	## year = identifies year for which ucr request was made

## Figure 1 ##
	deathsCount <- summaryBy(police_homicide + police_homicide_wp + 
		police_homicide_guardian + police_homicide_mpv +
		police_homicide_kbp + police_homicide_fatal ~ year, FUN=c(sum), data=data)

	deathsCount <- subset(deathsCount, year < 2017)
	deathsCountlong <- gather(deathsCount, source, deaths, -year)
	names(deathsCountlong)[names(deathsCountlong) == 'source'] <- 'Source'

	deathsCountlong$Source[deathsCountlong$Source=="police_homicide.sum" ] <- "FBI"
	deathsCountlong$Source[deathsCountlong$Source=="police_homicide_mpv.sum" ] <- "MPV"
	deathsCountlong$Source[deathsCountlong$Source=="police_homicide_fatal.sum" ] <- "Fatal"
	deathsCountlong$Source[deathsCountlong$Source=="police_homicide_kbp.sum" ] <- "KBP"
	deathsCountlong$Source[deathsCountlong$Source=="police_homicide_wp.sum" ] <- "WaPost"
	deathsCountlong$Source[deathsCountlong$Source=="police_homicide_guardian.sum" ] <- "Guardian"

	deathsCountlong$Source <- factor(deathsCountlong$Source,
		levels = c("FBI", "MPV", "Fatal", "KBP", "WaPost", "Guardian"))

## write the figure ##
		ggplot(deathsCountlong, aes(year, deaths, fill = Source)) +  
			geom_bar(position = "dodge", colour="black", size=0.25, stat="identity") +
			geom_text(aes(label=deaths), vjust = -0.25, size =1.5,fontface = "bold", position = position_dodge(width=0.9)) + 
			xlab("Year") + ylab("Police Killings")+ scale_fill_brewer(palette = "Greys", direction = -1) + 
	  		theme_bw()  
		ggsave("figure1.pdf", height=4, width=6, units="in", dpi=1000)
	
## Figure 2 ##
	deathsMap <- subset(data, year == 2015)
	deathsMap <- deathsMap[c("state_abb","police_homicide", "police_homicide_mpv", "police_homicide_fatal",
		"police_homicide_kbp", "police_homicide_wp", "police_homicide_guardian", "pop_100000")]
	deathsMap$police_homicide <- deathsMap$police_homicide/(deathsMap$pop_100000)
	deathsMap$police_homicide_mpv <- deathsMap$police_homicide_mpv/(deathsMap$pop_100000)
	deathsMap$police_homicide_fatal <- deathsMap$police_homicide_fatal/(deathsMap$pop_100000)
	deathsMap$police_homicide_kbp <- deathsMap$police_homicide_kbp/(deathsMap$pop_100000)
	deathsMap$police_homicide_wp <- deathsMap$police_homicide_wp/(deathsMap$pop_100000)
	deathsMap$police_homicide_guardian <- deathsMap$police_homicide_guardian/(deathsMap$pop_100000)

	deathsMap <- deathsMap[c("state_abb","police_homicide", "police_homicide_mpv", "police_homicide_fatal",
		"police_homicide_kbp", "police_homicide_wp", "police_homicide_guardian")]
	deathsMap <- gather(deathsMap, source, deaths, -state_abb)

	names(deathsMap)[names(deathsMap) == 'source'] <- 'Source'
	deathsMap$Source[deathsMap$Source=="police_homicide" ] <- "FBI"
	deathsMap$Source[deathsMap$Source=="police_homicide_mpv" ] <- "MPV"
	deathsMap$Source[deathsMap$Source=="police_homicide_fatal" ] <- "Fatal"
	deathsMap$Source[deathsMap$Source=="police_homicide_kbp" ] <- "KBP"
	deathsMap$Source[deathsMap$Source=="police_homicide_wp" ] <- "WaPost"
	deathsMap$Source[deathsMap$Source=="police_homicide_guardian" ] <- "Guardian"
	deathsMap$Source <-factor(deathsMap$Source,
		levels = c("FBI", "MPV", "Fatal", "KBP", "WaPost", "Guardian"))

## write the figure ##
		ggplot(deathsMap, aes(Source, deaths, fill = Source)) +
			#geom_bar(position = "dodge", colour="black", stat="identity") +
			geom_col(colour="black", size =0.25) +
		  theme_bw() +
			ylab("Police Killings/Population (per 100,000)") +
			facet_geo(~ state_abb, grid = "us_state_without_DC_grid2")  +
				scale_fill_brewer(palette = "Greys", direction = -1) +
			theme(strip.background = element_blank(),
		    axis.title.x=element_blank(),
    		axis.text.x=element_blank(),
    		axis.ticks.x=element_blank())
		ggsave("figure2.pdf", height=6, width=8, units="in", dpi=1000)
		

## Table 2 models ##
## washington post ##
	reWPdif <- lmer(wapost_fbi_dif/pop_100000 ~ std_squire_capacity + termEnact +
		std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov +
		demGov + std_black + (1|state), data = data)

## guardian ##
	reGdif <- lmer(guard_fbi_dif/pop_100000 ~ std_squire_capacity + termEnact +
		std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov +
		demGov + std_black + (1|state), data = data)

## mapping police violence ##
	reMPVdif <- lmer(mpv_fbi_dif/pop_100000 ~ std_squire_capacity + termEnact +
		std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov +
		demGov + std_black + (1|state), data = data)

## killed by police ##
	reKBPdif <- lmer(kbp_fbi_dif/pop_100000 ~ std_squire_capacity + termEnact +
		std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov +
		demGov + std_black + (1|state), data = data)

## fatal encounters ##
	reFATALdif <- lmer(fatal_fbi_dif/pop_100000 ~ std_squire_capacity + termEnact +
		std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov +
		demGov + std_black + (1|state), data = data)

## table 2 output ##
	stargazer(reWPdif, reGdif, reMPVdif, reKBPdif, reFATALdif, 
		title = "Unreported by FBI (Random Effects)",  
		dep.var.labels=c("WaPost", "Guardian", "MPV","KBP", "Fatal"), omit.stat=c("ser","f"))	

## table 3 models ##
## washington post ##
	lawWPdid <- lm(wapost_fbi_dif ~ invLaw + as.factor(state) + as.factor(year), data = data)

## guardian ##
	lawGdid <- lm(guard_fbi_dif ~ invLaw + as.factor(state) + as.factor(year), data = data)

## mapping police violence ##
	lawMPVdid <- lm(mpv_fbi_dif ~ invLaw + as.factor(state) + as.factor(year), data = data)

## killed by police ##
	lawKBPdid <- lm(kbp_fbi_dif ~ invLaw + as.factor(state) + as.factor(year), data = data)

## fatal encounters ##
	lawFATALdid <- lm(fatal_fbi_dif ~ invLaw + as.factor(state) + as.factor(year), data = data)

## table 3 output ##
	stargazer(lawWPdid, lawGdid, lawMPVdid, lawKBPdid, lawFATALdid,
		title = "Diff-in-Diff analysis of investigatory requirements",
		dep.var.labels=c("WaPost", "Guardian", "MPV","KBP", "Fatal"), omit.stat=c("ser","f"))

## now the appendix models / figures ##
## this is the Table A.1 analysis of legislative proposals ##
	proposals <- read.table('floydProposals6.24.csv', sep = ',', header = TRUE)
	proposals$state <- toupper(proposals$state)
	capacity <- unique(new[new$year == 2016, c('state', 'faOrd')])
	
	prop <- merge(proposals, capacity, by = c('state'), all.x = TRUE)
	
## MN has only split legislature, this recodes them such that ##
## we have Republican or otherwise ##
	prop$legPartyNew <- as.character(prop$legParty)
	prop$legPartyNew[prop$state == 'MINNESOTA'] <- 'dem'

	nb <- glm.nb(proposals ~ faOrd + legParty * govParty, data = prop)
	summary(nb)

	nbMnRecode <- glm.nb(proposals ~ faOrd + legPartyNew * govParty, data = prop)
	summary(nbMnRecode)

	nbMnOmit <- glm.nb(proposals ~ faOrd + legPartyNew * govParty, data = prop[prop$state != 'MINNESOTA', ])
	summary(nbMnOmit)
	
## now linear estimation of logged counts ##
	log <- lm(log(proposals + 1) ~ faOrd + legParty * govParty, data = prop)
	summary(log)

	logMnRecode <- lm(log(proposals + 1) ~ faOrd + legPartyNew * govParty, data = prop)
	summary(logMnRecode)

	logMnOmit <- lm(log(proposals + 1) ~ faOrd + legPartyNew * govParty, data = prop[prop$state != 'MINNESOTA', ])
	summary(logMnOmit)

## output Table A.1 ##
	stargazer(nbMnRecode, nbMnOmit, logMnRecode, logMnOmit)
	
## Figure A.1 plot of DV ##
	all <- rbind(old, new)
	all <- all[all$year != 1977, ]
	all <- summaryBy(ucr ~ state + year, keep.names = TRUE, data = all)
	means <- summaryBy(ucr ~ state, FUN = mean, data = all, keep.names = TRUE)
	colnames(means)[2] <- 'means'
	hold <- merge(all, means, by = 'state', all.x = TRUE)
	hold <- hold[order(hold$means), ]
	newOrder <- with(hold, reorder(state, ucr, median))
	
	pdf('compliance.pdf', width = 10)
	boxplot(hold$ucr ~ newOrder, 
		ylim = c(-0.3, 1),
		axes = FALSE,
		col = rgb(0, 0, 0, 0.2),
		border = rgb(0, 0, 0, 0.5),
		main = 'UCR Compliance by State',
		xlab = 'States (Median Ranked)',
		ylab = 'UCR Compliance')
	axis(2, at = c(0, 1))
	for(i in 1:length(unique(hold$state))){
		j <- unique(hold$state)[i]
		
		text(i, -0.03, j, 
			srt = 90, adj = 1, cex = 0.75, col = rgb(0, 0, 1, 1))
	}
	dev.off()

## Figure A.2, plot of IV ##
	all <- rbind(old, new)
	all <- summaryBy(faOrd ~ state + year, keep.names = TRUE, data = all)
	means <- summaryBy(faOrd ~ state, FUN = mean, data = all, keep.names = TRUE)
	colnames(means)[2] <- 'means'
	hold <- merge(all, means, by = 'state', all.x = TRUE)
	hold <- hold[order(hold$means), ]
	newOrder <- with(hold, reorder(state, faOrd, median))
	
	pdf('capacity.pdf', width = 10)
	boxplot(hold$faOrd ~ newOrder, 
		ylim = c(-5, 4),
		axes = FALSE,
		col = rgb(0, 0, 0, 0.2),
		border = rgb(0, 0, 0, 0.5),
		main = 'Legislative Capacity by State',
		xlab = 'States (Median Ranked)',
		ylab = 'Legislative Capacity')
	axis(2)
	for(i in 1:length(unique(hold$state))){
		j <- unique(hold$state)[i]
		
		text(i, min(hold$faOrd[hold$state == j])-0.4, j, 
			srt = 90, adj = 1, cex = 0.75, col = rgb(0, 0, 1, 1))
	}
	dev.off()

## Table A.3 plot of term limit imposition ##
	all <- rbind(old, new)
	all <- summaryBy(cbind(termEnact, ucr) ~ state + year, keep.names = TRUE, data = all)
	pdf('enact.pdf', width = 10, height = 8)
	panelView(ucr ~ termEnact, data = all, index = c('state', 'year'), pre.post = TRUE,
		main = 'Term Limit Enactment by State Year', xlab = 'Year', ylab = 'State', 
		background = 'white', axis.lab.gap = c(2, 0))
	dev.off()
	
## Table A.2 is descriptives ##
	summary(old)
	summary(new)

## Table A.3 is the full results of the main compliance models ##
	stargazer(squireOld, termOld, fullOld, squireNew, termNew, fullNew)

## Table A.4 compares our squire measure to the VAM measure ##
	sumNew <- summaryBy(cbind(vamDv, faOrd) ~ state + year, data = new)
	modNew <- lm(vamDv.mean ~ faOrd.mean, data = sumNew)
	summary(modNew)
	sumOld <- summaryBy(cbind(vamDv, faOrd) ~ state + year, data = old)
	modOld <- lm(vamDv.mean ~ faOrd.mean, data = sumOld)
	summary(modOld)
	stargazer(modOld, modNew)

## Table A.5 compares DiD results between squire and VAM measure ##
## estimate the VAM models ##
	vamOld <- lm(ucr ~ vamDv + 
		as.factor(state) + as.factor(year), data = old)
	summary(vamOld)

	vamOldC <- lm(ucr ~ (vamDv + termEnact) * election + 
		type + conSheriff + scale(employees) + scale(population) + 
		termEnact + demGov + repGov + black + 
		as.factor(state) + as.factor(year)*(demLeg + repLeg), data = old)
	summary(vamOldC)

	vamNew <- lm(ucr ~ vamDv + 
		as.factor(state) + as.factor(year), data = new)

	vamNewC <- lm(ucr ~ (vamDv + termEnact) * election + 
		type + conSheriff + scale(employees) + scale(population) + 
		termEnact + demGov + repGov + black + 
		as.factor(state) + as.factor(year)*(demLeg + repLeg), data = new)
	summary(vamNewC)
	stargazer(squireOld, vamOld, fullOld, vamOldC, squireNew, vamNew, fullNew, vamNewC)

## Table A.6 breaks models into city/county/state levels ##
	oldCity <- lm(ucr ~ faOrd + termEnact +
		as.factor(state) + as.factor(year), data = old[old$type == 'City', ])
	summary(oldCity)

	oldCounty <- lm(ucr ~ faOrd + termEnact +
		as.factor(state) + as.factor(year), data = old[old$type == 'County', ])
	summary(oldCounty)

	oldState <- lm(ucr ~ faOrd + termEnact +
		as.factor(state) + as.factor(year), data = old[old$type == 'State Police', ])
	summary(oldState)

	newCity <- lm(ucr ~ faOrd + termEnact +
		as.factor(state) + as.factor(year), data = new[new$type == 'City', ])
	summary(newCity)

	newCounty <- lm(ucr ~ faOrd + termEnact +
		as.factor(state) + as.factor(year), data = new[new$type == 'County', ])
	summary(newCounty)

	newState <- lm(ucr ~ faOrd + termEnact +
		as.factor(state) + as.factor(year), data = new[new$type == 'State Police', ])
	summary(newState)

	stargazer(oldCity, oldCounty, oldState, newCity, newCounty, newState)

## Table A.7 allows effects to vary by party control of legislature ##
	squireOldC <- lm(ucr ~ (faOrd + termEnact) * (election + demLeg + repLeg) + 
		conSheriff + type + scale(employees) + scale(population) + black + 
		termEnact + demGov + repGov + 
		as.factor(state) + as.factor(year)*(demLeg + repLeg), data = old)
	summary(squireOldC)
repLeg + demLeg 
		conSheriff + type + scale(employees) + scale(population) + black + 
		termEnact + demGov + repGov + 
		as.factor(state) + as.factor(year)*(demLeg + repLeg), data = new)
	summary(squireNewC)

	stargazer(squireOld, termOld, squireOldC, squireNew, termNew, squireNewC)

## the next set of results focuses on the underreporting analyses ##
## Table A.8 seperate outcomes ##
	FBIall <- lmer(police_homicide/pop_100000 ~  std_squire_capacity  + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov +  (1|state), data = data)
	WPall <- lmer(police_homicide_wp/pop_100000 ~ std_squire_capacity  + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov +  (1|state), data = data)
	Gall <- lmer(police_homicide_guardian/pop_100000 ~ std_squire_capacity  + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov +  (1|state), data = data)
	MPVall <- lmer(police_homicide_mpv/pop_100000 ~  std_squire_capacity  + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov +  (1|state), data = data)
	KBPall <- lmer(police_homicide_kbp/pop_100000 ~  std_squire_capacity  + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov +  (1|state), data = data)
	FATALall <- lmer(police_homicide_fatal/pop_100000 ~ std_squire_capacity  + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov +  (1|state), data = data)

	stargazer(FBIall, WPall, Gall, MPVall, KBPall, FATALall, title = "Reported Deaths (OLS) -- All",
		dep.var.labels=c("FBI","WaPost", "Guardian", "MPV","KBP", "Fatal"), omit.stat=c("ser","f"))	

## time FE Table A.9 ##
	data$year_factor <- as.factor(data$year)
	data <- within(data, year_factor <- relevel(year_factor, ref = "2016"))
	WPdif_timeFE <- lmer(wapost_fbi_dif/pop_100000 ~ std_squire_capacity  + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov + 
		as.factor(year_factor) + (1|	state), data = data)
	Gdif_timeFE <- lmer(guard_fbi_dif/pop_100000 ~ std_squire_capacity  + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov + 
		as.factor(year_factor) + (1| state), data = data)
	MPVdif_timeFE <- lmer(mpv_fbi_dif/pop_100000 ~  std_squire_capacity  + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov + 
		as.factor(year_factor) + (1| state), data = data)
	KBPdif_timeFE <- lmer(kbp_fbi_dif/pop_100000 ~  std_squire_capacity  + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov + 
		as.factor(year_factor) + (1| state), data = data)
	FATALdif_timeFE <- lmer(fatal_fbi_dif/pop_100000 ~ std_squire_capacity  + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov + 
		as.factor(year_factor) + (1| state), data = data)

	stargazer(WPdif_timeFE, Gdif_timeFE, MPVdif_timeFE, KBPdif_timeFE, FATALdif_timeFE, title =  "Unreported by FBI (Random Effects) -- w/ Time Fixed Effects",
		dep.var.labels=c("WaPost", "Guardian", "MPV","KBP", "Fatal"), omit.stat=c("ser","f"))	

## now year-by-year models ##
## 2013 Table A.10 ##
	MPVdif2013 <- lm(mpv_fbi_dif/pop_100000 ~  std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov, data = subset(data, year == 2013))
	FATALdif2013 <- lm(fatal_fbi_dif/pop_100000 ~  std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov, data = subset(data, year == 2013))
	
	stargazer(MPVdif2013, FATALdif2013, title =  "Unreported by FBI - 2013",  dep.var.labels=c( "MPV", "Fatal"), omit.stat=c("ser","f"))	

## 2014 A.11 ##
	KBPdif2014 <- lm(kbp_fbi_dif/pop_100000 ~ std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov, data = subset(data, year == 2014))
	MPVdif2014 <- lm(mpv_fbi_dif/pop_100000 ~ std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov, data = subset(data, year == 2014))
	FATALdif2014 <- lm(fatal_fbi_dif/pop_100000 ~ std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov, data = subset(data, year == 2014))

	stargazer(KBPdif2014,MPVdif2014, FATALdif2014, title =  "Unreported by FBI - 2014",  dep.var.labels=c("KBP", "MPV", "Fatal"), omit.stat=c("ser","f"))	

## 2015 A.12 ##
	WPdif2015 <- lm(wapost_fbi_dif/pop_100000 ~ std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov, data = subset(data, year == 2015))
	Gdif2015 <- lm(guard_fbi_dif/pop_100000 ~ std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov, data = subset(data, year == 2015))
	KBPdif2015 <- lm(kbp_fbi_dif/pop_100000 ~ std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov, data = subset(data, year == 2015))
	MPVdif2015 <- lm(mpv_fbi_dif/pop_100000 ~ std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov, data = subset(data, year == 2015))
	FATALdif2015 <- lm(fatal_fbi_dif/pop_100000 ~ std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov, data = subset(data, year == 2015))

	stargazer(WPdif2015,Gdif2015,KBPdif2015,MPVdif2015, FATALdif2015, title =  "Unreported by FBI - 2015",
		dep.var.labels=c("WaPost", "Guardian", "MPV","KBP", "Fatal"), omit.stat=c("ser","f"))	

## 2016 A.13 ##
	WPdif2016 <- lm(wapost_fbi_dif/pop_100000 ~ std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov, data = subset(data, year == 2016))
	Gdif2016 <- lm(guard_fbi_dif/pop_100000 ~ std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov, data = subset(data, year == 2016))
	KBPdif2016 <- lm(kbp_fbi_dif/pop_100000 ~ std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov, data = subset(data, year == 2016))
	MPVdif2016 <- lm(mpv_fbi_dif/pop_100000 ~ std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov, data = subset(data, year == 2016))
	FATALdif2016 <- lm(fatal_fbi_dif/pop_100000 ~ std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov, data = subset(data, year == 2016))

	stargazer(WPdif2016,Gdif2016,KBPdif2016,MPVdif2016, FATALdif2016, title =  "Unreported by FBI - 2016",
		dep.var.labels=c("WaPost", "Guardian", "MPV","KBP", "Fatal"), omit.stat=c("ser","f"))	

## count distribution models ##
## censor the data from below ##
	data$wapost_fbi_dif_censor <- ifelse(data$wapost_fbi_dif<0, 0, data$wapost_fbi_dif)
	data$guard_fbi_dif_censor <- ifelse(data$guard_fbi_dif<0, 0, data$guard_fbi_dif)
	data$mpv_fbi_dif_censor <- ifelse(data$mpv_fbi_dif<0, 0, data$mpv_fbi_dif)
	data$kbp_fbi_dif_censor <- ifelse(data$kbp_fbi_dif<0, 0, data$kbp_fbi_dif)
	data$fatal_fbi_dif_censor <- ifelse(data$fatal_fbi_dif<0, 0, data$fatal_fbi_dif)

## poisson Table A.14 ##
	pWPdif <- glm(wapost_fbi_dif_censor ~ std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov + 
		log(pop_100000), family = "poisson", data = data)
	pGdif <- glm(guard_fbi_dif_censor ~ std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov + 
		log(pop_100000), family = "poisson", data = data)
	pMPVdif <- glm(mpv_fbi_dif_censor ~  std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov + 
		log(pop_100000) , family = "poisson", data = data)
	pKBPdif <- glm(kbp_fbi_dif_censor ~  std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov + 
		log(pop_100000) ,  family = "poisson", data = data)
	pFATALdif <- glm(fatal_fbi_dif_censor ~ std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov + 
		log(pop_100000),  family = "poisson", data = data)
	
	stargazer(pWPdif, pGdif, pMPVdif, pKBPdif, pFATALdif, title =  "Unreported by FBI (Poisson)", 
		dep.var.labels=c("WaPost", "Guardian", "MPV","KBP", "Fatal"), omit.stat=c("ser","f"))	

## negative binomial Table A.15 ##
	nbWPdif <- glm.nb(wapost_fbi_dif_censor ~ std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov + 
		log(pop_100000), data = data)
	nbGdif <- glm.nb(guard_fbi_dif_censor ~ std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov + 
		log(pop_100000), data = data)
	nbMPVdif <- glm.nb(mpv_fbi_dif_censor ~ std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov + 
		log(pop_100000), data = data)
	nbKBPdif <- glm.nb(kbp_fbi_dif_censor ~ std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov + 
		log(pop_100000), data = data)
	nbFATALdif <- glm.nb(fatal_fbi_dif_censor ~  std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov + 
		log(pop_100000), data = data)

	stargazer(nbWPdif, nbGdif, nbMPVdif, nbKBPdif, nbFATALdif, title =  "Unreported by FBI (NB)", 
		dep.var.labels=c("WaPost", "Guardian", "MPV","KBP", "Fatal"), omit.stat=c("ser","f"))	

## now the investigation law content ##
## balance picture --- Figure A.4 ##
pdf('balance.pdf')
	panelView(mpv_fbi_dif ~ invLaw, data = data, index = c('state', 'year'), pre.post = TRUE,
		main = 'Enactment of Investigation Laws', xlab = 'Year', ylab = 'State', background = 0)
dev.off()

## Table A.16 is just the same as the main text ##
	stargazer(lawWPdid, lawGdid, lawMPVdid, lawKBPdid, lawFATALdid,
		title = "Diff-in-Diff analysis of investigatory requirements",
		dep.var.labels=c("WaPost", "Guardian", "MPV","KBP", "Fatal"), omit.stat=c("ser","f"))

## Table A.17 adds control variables ##

	lawWPdidC <- lm(wapost_fbi_dif ~ invLaw +
		pop_100000 + std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov + 
		as.factor(state) + as.factor(year), data = data)
	lawGdidC <- lm(guard_fbi_dif ~ invLaw +
		pop_100000 + std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov + 
		as.factor(state) + as.factor(year), data = data)
	lawMPVdidC <- lm(mpv_fbi_dif ~ invLaw +
		pop_100000 + std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov + 
		as.factor(state) + as.factor(year), data = data)
	lawKBPdidC <- lm(kbp_fbi_dif ~ invLaw +
		pop_100000 + std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov + 
		as.factor(state) + as.factor(year), data = data)
	lawFATALdidC <- lm(fatal_fbi_dif ~ invLaw +
		pop_100000 + std_squire_capacity + std_med_income + std_employees + std_police_ideal + repLeg + demLeg + repGov + demGov + 
		as.factor(state) + as.factor(year), data = data)

	stargazer(lawWPdidC, lawGdidC, lawMPVdidC, lawKBPdidC, lawFATALdidC, title =  "Diff-in-Diff analysis of investigatory requirements",
		dep.var.labels=c("WaPost", "Guardian", "MPV","KBP", "Fatal"), omit.stat=c("ser","f"))

## the LEMAS replication given in Table A.18 ##
	lemas <- read.table('lemas2016.txt', sep = '|', header = TRUE)
	hold <- new[new$year == 2016, ]
	hold <- hold[, c('state', 'faOrd', 'termEnact')]
	hold <- unique(hold)
	lemas$state <- toupper(state.name[match(lemas$state, state.abb)])
	lemas <- merge(hold, lemas, by = 'state')
	summary(lemas)

	base <- lm(completed ~ faOrd, data = lemas)
	term <- lm(completed ~ termEnact, data = lemas)
	both <- lm(completed ~ faOrd + termEnact, data = lemas)
	stargazer(base, term, both)

## are constitutional sheriffs really different? ##
## wyoming constitutional change ##
## Figure A.5 ##

	table(old$year[old$state == 'WYOMING' & old$type == 'County'],
		old$ucr[old$state == 'WYOMING' & old$type == 'County'])
	table(new$year[new$state == 'WYOMING' & new$type == 'County'],
		new$ucr[new$state == 'WYOMING' & new$type == 'County'])
	years <- c(old$year[old$state == 'WYOMING' & old$type == 'County'],
		new$year[new$state == 'WYOMING' & new$type == 'County'])
	ucr <- c(old$ucr[old$state == 'WYOMING' & old$type == 'County'],
		new$ucr[new$state == 'WYOMING' & new$type == 'County'])
	wy <- data.frame(years, ucr)
	wy <- summaryBy(ucr ~ years, FUN = mean, data = wy, keep.names = TRUE)

	pdf('picWy.pdf')
	plot(wy$years, wy$ucr, type = 'n',
		xlab = 'Year',
		ylim = c(0.75, 1),
		ylab = 'Compliance Rate',
		main = 'Wyoming: Redaction of Sheriff from Constitution',
		axes = FALSE)
	axis(1)
	axis(2)
	points(wy$years[wy$years < 1991], wy$ucr[wy$years < 1991], col = rgb(0, 0, 1, 0.7), pch = 16)
	lines(lowess(wy$years[wy$years < 1991], wy$ucr[wy$years < 1991]), col = rgb(0, 0, 1, 0.7), lwd = 2)
	
	points(wy$years[wy$years >= 1991], wy$ucr[wy$years >= 1991], col = rgb(1, 0, 0, 0.7), pch = 17)
	lines(lowess(wy$years[wy$years >= 1991], wy$ucr[wy$years >= 1991]), col = rgb(1, 0, 0, 0.7), lwd = 2)
	
	lines(c(1990.5, 1990.5), c(0, 1))
	text(2000, 0.85, 'Constitutional reform')
	dev.off()

## Table A.20 models ##
	years <- c(old$year[old$state == 'WYOMING'],
		new$year[new$state == 'WYOMING'])
	ucr <- c(old$ucr[old$state == 'WYOMING'],
		new$ucr[new$state == 'WYOMING'])
	type <- c(old$type[old$state == 'WYOMING'],
		new$type[new$state == 'WYOMING'])
	name <- c(old$ori[old$state == 'WYOMING'],
		new$ori[new$state == 'WYOMING'])
	wy <- data.frame(years, ucr, type, name)

	mod <- lm(ucr ~ (years > 1990)*(type == 'County') +as.factor(years), data = wy)
	summary(mod)	
	stargazer(mod)

## what about within-unit differences in alabama? ##
## geneva, lamar, pickens, and st. clare counties have funding gurantees ##
## let's have a look ##
## those are ori: AL0340000, AL0400000, AL0540000, AL0580000 ##

	subOld <- old[old$state == 'ALABAMA' & old$type == 'County', ]
	subOld$fund <- 0
	subOld$fund[subOld$ori == 'AL0340000' | subOld$ori == 'AL0400000' |
		subOld$ori == 'AL0540000' | subOld$ori == 'AL0580000'] <- 1
	subNew <- new[new$state == 'ALABAMA' & new$type == 'County', ]
	subNew$fund <- 0
	subNew$fund[subNew$ori == 'AL0340000' | subNew$ori == 'AL0400000' |
		subNew$ori == 'AL0540000' | subNew$ori == 'AL0580000'] <- 1

	years <- c(subOld$year, subNew$year)
	ucr <- c(subOld$ucr, subNew$ucr)
	fund <- c(subOld$fund, subNew$fund)
	al <- data.frame(years, ucr, fund)
	t.test(al$ucr[al$fund == 0], al$ucr[al$fund == 1])
	al <- summaryBy(ucr ~ years + fund, FUN = mean, data = al, keep.names = TRUE)

	table(subOld$ucr, subOld$fund)
	table(subNew$ucr, subNew$fund)

## this is Figure A.6 ##
	pdf('picAl.pdf')
	plot(al$years, al$ucr, type = 'n',
		xlab = 'Year',
		ylim = c(0.25, 1),
		ylab = 'Compliance Rate',
		main = 'Alabama: Constitutional fund mandates',
		axes = FALSE)
	axis(1)
	axis(2)
	points(al$years[al$fund == 0], al$ucr[al$fund == 0], col = rgb(0, 0, 1, 0.7), pch = 16)
	lines(lowess(al$years[al$fund == 0], al$ucr[al$fund == 0]), col = rgb(0, 0, 1, 0.7), lwd = 2)
	
	points(al$years[al$fund == 1], al$ucr[al$fund == 1], col = rgb(1, 0, 0, 0.7), pch = 17)
	lines(lowess(al$years[al$fund == 1], al$ucr[al$fund == 1]), col = rgb(1, 0, 0, 0.7), lwd = 2)
	
	legend('bottomright', legend = c('Fund Mandate', 'No Fund Mandate'), 
		col = c(rgb(1, 0, 0, 0.7), rgb(0, 0, 1, 0.7)), pch = c(17, 16))
	dev.off()
	
